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A general algorithm toward the solution of the fermion sign 
problem in finite-temperature quantum Monte Carlo simula- 
tions has been formulated for discretized fermion path in- 
tegrals with nearest-neighbor interactions in the Trotter di- 
rection. This multilevel approach systematically implements 
a simple blocking strategy in a recursive manner to synthe- 
size the sign cancellations among different fermionic paths 
throughout the whole configuration space. The practical use- 
fulness of the method is demonstrated for interacting electrons 
in a quantum dot. 

PACS numbers: 02.70.Lq, 05.30.Fk, 73.20.Dx 



The quantum Monte Carlo (QMC) technique is one of 
the most powerful methods for the simulation of many- 
fermion systems. It is based on a path integral for- 
mulation of the fermion propagator and is one of the 
very few methods capable of delivering exact results for 
strongly correlated systems. Despite its potentials, ap- 
plications of QMC have been severely handicapped by 
the notorious "fermion sign problem" As a con- 

sequence of exchange, fermionic density matrix elements 
are not positive-definite. The sign cancellations arising 
from sampling fermion paths then manifest themselves as 
a small signal-to-noise ratio that vanishes exponentially 
with either the system size or with decreasing tempera- 
ture. Besides variational or approximate treatments such 
as the fixed-node approximation j^] , the sign problem has 
remained unsolved. 

In this Letter, we propose a simple and intuitive ap- 
proach toward the general solution of the fermion sign 
problem. Our algorithm represents the systematic imple- 
mentation of a blocking strategy [0 . The idea behind the 
blocking strategy is that by sampling groups of states, the 
sign problem can always be reduced compared to sam- 
pling single states. By suitably bunching states together 
into blocks, the sign cancellations among states within 
the same block can be accounted for non-stochastically. 
It can then be shown [Q that any such blocking will al- 
ways reduce the sign problem — no blocking will ever 
make the sign problem more severe 0j . Any real progress 
on the sign problem will require an accurate treatment 
of the sign cancellations within suitably chosen subunits 
(blocks) of state space. 

A systematic improvement of the sign problem can be 



achieved by formulating the blocking strategy in a recur- 
sive bottom-to-top fashion. Blocks of different sizes are 
defined on several levels, and after taking care of the sign 
cancellations within all blocks on a given (finer) level, 
the resulting sign problem can be transferred to the next 
(coarser) level. By doing this recursively, the sign prob- 
lem on all the coarser levels can be handled in the same 
manner. It is then possible to proceed without numerical 
instabilities from the bottom up to the top level, where 
the last remaining cancellations pose no serious challenge. 

In many ways, the algorithm we are proposing is re- 
lated to the renormalization group approach. But instead 
of integrating out information on fine levels, the sign can- 
cellations are "synthesized" within a given level and sub- 
sequently their effects are transferred to coarser levels. 
Our approach is actually closer in spirit to the multi-grid 
algorithm [|| . The method of Rcf . jl| can be understood 
as a uni-level scheme which provides a partial solution 
to the sign problem. In contrast, given sufficient com- 
puter memory, the algorithm proposed here can provide 
a complete solution. Below we describe this multilevel 
blocking (MLB) algorithm and apply it to the simulation 
of correlated electrons in a quantum dot. 

We consider a many-fermion system whose state is de- 
scribed by a set of quantum numbers r denoting, e.g., 
the positions and spins of all particles. For simplicity, we 
focus on calculating the expectation value of a diagonal 
operator 0, 



T ir A ( r )p( r , r ) 



(1) 



where ^ r represents either a summation for the case of a 
discrete system or an integration for a continuous system. 
Imaginary time is then discretized into P slices of length 
e = (3/P, where (3 — \/ksT and we require P = 2 L . 
Inserting complete sets at each slice rn = 1, . . . , P, and 
denoting the corresponding configuration on slice m by 
r m , the diagonal elements of the density matrix at r ~ 
Tp read 



p(P,P)= ]T (l,2) (2,3)o---(P,i; 



o ■ 



(2) 



As a shorthand notation, we use the slice index m for 
the quantum numbers r m . This equation also defines 
the level- bonds, which are simply the short-time prop- 
agators, 
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(m,m + l)o = (r m+ i\e e \r m ) 



(3) 



This formulation of the problem excludes effective ac- 
tions such as those arising from an integration over the 
fcrmions via the Hubbard-Stratonovich transformation 
, since they generally lead to long-ranged imaginary- 
time interactions. 
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FIG. 1. Levels for L = 2 (P = 4). Imaginary time flows 
along the circle (solid curve), and the slices m = 1, 2, 3, 4 are 
distributed among the three levels: The finest level I — 
contains m = 1,3, level I = 1 contains m = 2, and £ = 2 
contains m = 4. Levels bonds are indicated by dashed and 
dotted lines. 

To describe the MLB strategy, we need to specify the 
different levels < £ < L, where L defines the Trotter 
number P = 2 L . Each slice m belongs to a unique level 
£, such that m = (2j + l)2 l and j is a nonnegative inte- 
ger. For instance, the slices m = 1, 3, 5, • ■ • , P — 1 belong 
to t = 0, m = 2, 6, 10, • • • , P - 2 belong to £ = 1, etc., 
such that there are A^ = 2 i_£_1 (but A/i, = 1) different 
slices on level £, see Figure [i] An elementary blocking is 
achieved by grouping together configurations that differ 
only at slice m, so only r m varies in that block while all 
r m i^ m remain fixed. Sampling on level £ therefore ex- 
tends over configurations {r m } living on the Mi different 
slices. In the MLB scheme, we move recursively from the 
finest (£ = 0) up to the coarsest level (£ = L), and the 
measurement of the diagonal operator is done only at the 
top level using the configuration rp. 

We now describe a practical implementation of the 
MLB scheme. A Monte Carlo sweep starts by chang- 
ing only configurations associated with the slices on level 
£ = according to the weight 



P = |(l,2)o(2,3)o---(P,l)o| 



(4) 



generating a MC trajectory containing K samples for 
each slice on level £ = 0. These NqK samples are stored 
and they are used to generate additional coarser interac- 
tions among the higher-level slices, 



(m, m + 2)i = (sgn[(m, m + l)o(m + 1, m + 2) ]}p [ m+ i] 

m+l 

x (m + l,?7i + 2) ] , 

where the summation X)m+i extends over the MqK sam- 
ples. As will be discussed in detail later on, for a com- 
plete solution of the sign problem, the sample number K 
should be chosen as large as possible. The level-1 bonds 
(||) contain crucial information about the sign cancella- 
tions on the previous level £ = 0. Using these bonds, the 
density matrix (^|) is rewritten as 



P(P,P) 



E 

1,2,...,P-1 



|(l,2) (2,3)o-"(P,l)o| 



x (2, 4)i • • • (P — 2, P)i(P, 2)i . 



(0) 



Comparing this to Eq. (||), we notice that the entire sign 
problem has been transferred to the next coarser level by 
using the level-1 bonds. 

In the next step, the sampling is carried out on level 
I = 1 in order to generate the next-level bonds, i.e., only 
slices m = 2, 6, . . . , P — 2 are updated, using the weight 
VqPi with 



V 1 = |(2,4)i(4,6)i--.(P,2)i| 



(7) 



Moving the level-1 configurations modifies the level-0 
bonds, which in turn requires that the level-1 bonds be 
updated. A direct re-calculation of these bonds accord- 
ing to Eq. (||) would be too costly. Instead, we use the 
stored configurations on level £ = to perform an impor- 
tance sampling of the new level-1 bonds. Under the test 
move m — > m! (i.e., r m — > r' m ) on level 1=1, the bond 
(BO can be obtained from 



(m', 



E(m ,m+l)o(rn+i??7i'+2)o 
m+l |(rn 1 m+l)o(m+l 1 m+2)o| 
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(m',m+l)o(ro+l,TO+2) | 
m+l \(m,m+l)a(m+l,m+2) a \ 
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where Em+i runs over the previously stored MC config- 
urations r m+1 . Note that for small values of AT, Eq. (||) is 
only approximative, and thus a sufficiently large value of 
K should be chosen. With the aid of Eq. (JsJ) , we obtain 
the updated level-1 bonds with only moderate compu- 
tational effort. Generating a sequence of K samples for 
each slice on level 1=1, and storing these M\K samples, 
we then calculate the level-2 bonds in analogy to Eq. (|^), 

(m,m + 4) 2 = (sgn [(777,777 + 2)i(m + 2, m + 4)i]) Pl -p , 

(9) 

and iterate the process up to the top level £ = L using 
the obvious recursive generalization of Eqs. (JsJ) and (|^) 
to define level-£ bonds. 

Thereby the diagonal elements of the density matrix 
are obtained as 
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p{P,P)= Y, l(l,2)o(2,3)o--.(P,l)o| (10) 

1,2 P-l 

x 1(2,4)! •••(P-2,P) 1 (P,2) 1 | 
•■■|(P/2,P) L _ 1 (P,P/2) L _ 1 |(P,P) L . 

By virtue of this algorithm, the sign problem is trans- 
ferred step by step up to the coarsest level. The expec- 
tation value (GJ) can thus be computed from 



(A) = 



(A(P) sgn(P,P) L ) P 
(aga(P,P) L )v 



(11) 



The manifestly positive definite MC weight V used for 
the averaging in Eq. ( pi] ) can be read off from Eq. (|l0|), 



V = |(l,2) (2,3)o---(P,l)o| 
x 1(2,4)! •••(P-2,P) 1 (P,2) 1 | 
•••|(P/2,P) L _ 1 (P,P/2) i _ 1 ||(P,P) i | 



The denominator in Eq. ( pL 1| ) gives the average sign 
and indicates to what extent the sign problem has been 
solved. Under a naive application of the QMC technique, 
the average sign decays exponentially with (3 and is typi- 
cally close to zero. This causes the numerical instabilities 
associated with the sign problem, i.e., to obtain statis- 
tically relevant results requires exponentially long CPU 
times. With the MLB algorithm, however, the sign prob- 
lem can be completely eliminated. The average sign re- 
mains close to unity for any /3, with a CPU time require- 
ment that increases only linearly. The price to pay for 
the stability of the algorithm is the increased memory 
requirement associated with having to store the sampled 
configurations on the fine levels, which scales at worst 
quadratically in K. The example below demonstrates 
that modest memory requirements are sufficient to treat 
rather complex problems. 

Next we address questions concerning the exactness of 
the MLB approach for a finite sample number K. Clearly, 
K needs to be sufficiently large to produce a reliable es- 
timate for the levels bonds. If these bonds could be cal- 
culated exactly (corresponding to the limit K — > oo), the 
manipulations leading to Eq. ( |Io| ) yield the exact result. 
Hence for large enough K . the MLB algorithm must (i) 
become exact and (ii) completely solve the sign problem. 
However, since the level-^ bonds can only be computed 
for finite K 1 the weight function V amounts to using a 
noisy estimator, which in turn can introduce bias into the 
algorithm B. In principle, this problem could be avoided 
by using a linear acceptance criterion Q instead of the 
algorithmically simpler Metropolis choice M. But even 
with the Metropolis choice (which we used in the example 
below) , the bias can be made arbitrarily small by increas- 
ing K. Therefore, with sufficient computer memory, the 
MLB approach can be made to give numerically exact re- 
sults. One might then worry about the actual value of K 
required to obtain stable and exact results. If this value 



were to scale exponentially with j3 and/or system size, 
the sign problem would be present in disguise again. Al- 
though we do not have a rigorous non-exponential bound 
on K , our experience with the MLB algorithm indicates 
that this scaling is at worst algebraic. 

We now illustrate the usefulness of the method for 
N interacting electrons confined in a quantum dot [p|. 
Quantum dots are two-dimensional artificial atoms fab- 
ricated by means of suitable gates in semiconductor het- 
erostructures. For simplicity, we consider only spinlcss 
electrons, zero magnetic field, and a parabolic confine- 
ment potential. We employ a symmetric Trotter breakup 
for H = H 1 +H 2 , 



(12) Hi 



N ( v 2 

y — 




N 



h 2 = j:- 



i<3 



(13) 



where the positions (momenta) of the N electrons are 
Xj (pj), the dielectric constant is k, and m* is the ef- 
fective mass. The short-time propagator (^) under Hi 
is obtained by antisymmetrizing a product of N har- 
monic oscillator propagators, leading to a fermion de- 
terminant. Since the determinant can change sign, con- 
ventional QMC simulations run into the sign problem. 
The MC updating then employs randomly chosen single 
particle moves on the momentary level £. This suffices for 
an efficient and ergodic sampling of configuration space. 
Here we present results for the energy En = (H). Since 
the direct evaluation of the kinetic energy would involve a 
nonlocal operator, we have exploited the quantum virial 
theorem |jic|] in order to sample En. Simulations were 
done at Ti(3ujq = 6 for two different interaction strengths, 
lo/a = 2 and 8. Here lo = (h/m*^) 1 / 2 is a confinement 
lengthscale, and lo/a = e 2 / \%kIqujq) , with a being the 
effective Bohr radius. Trotter convergence was achieved 
at L = 3 (4) for l /a = 2 (8V The N = 2 exact di- 
agonalization results of Rcf. |11| have been accurately 



reproduced, which also serves as an independent check 
for our code. The simulations have been carried out on 
an IBM FJSC6000/590 workstation. 

TABLE I. MLB results for N = 8 and l /a = 2. N a is 
the number of samples (in 10 4 ), tepv the total CPU time (in 
hours), MB the required memory (in mega- bytes), and (sgn) 
the average sign. Bracketed numbers are error estimates. 
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MB 


(sgn) 


E N /hu 


1 


120 


95 


2 


0.02 


48.6(3) 


100 


7 


33 


14 


0.48 


48.43(8) 


200 


9 


83 


30 


0.63 


48.55(7) 


400 


8 


174 


64 


0.73 


48.53(9) 


600 


10 


308 


96 


0.77 


48.54(8) 


800 


9 


429 


150 


0.81 


48.59(8) 
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To elucidate how the MLB algorithm works in prac- 
tice, in Table | we compare simulation results for TV = 8 
electrons at various values of K. Compared to the naive 
approach (K = 1), using a moderate K = 200 already 
increases the average sign from 0.02 to 0.63, making it 
possible to obtain more accurate results from much fewer 
samples. The data in Table Q also confirms that the bias 
can be systematically eliminated by increasing K , so that 
the energy found at K > 200 essentially represents the 
exact result (within error bars). As expected, the CPU 
time per sample scales only linearly with K, where the 
memory requirements grow at most quadratically with 
K, 

Results for En with N < 8 are shown in Figure [|. 
For N < 5, the fixed- node QMC and Jastrow wavefunc- 
tion calculations of Ref. |l2| are in fairly good agreement 
with our exact results. However, for larger N, there are 
deviations, with the correct energies significantly lower 
than the values reported in Ref. [([3, which represent 
variational upper bounds. Notably, there are no obvious 
cusps or breaks in the iV-dependence of the energy. Such 
features would hint at the existence of magic numbers 
for which the artificial atom is exceptionally stable. Our 
data in Fig. |^ suggests that an explanation of the experi- 
mentally observed magic numbers || has to involve spin 
or magnetic field effects. Remarkably, the absence of pro- 
nounced cusps in Em/N for strong correlations (lo/a = 8) 
is in accordance with a purely classical analysis [Ofl . 




ized by local imaginary-time interactions. Given suffi- 
cient computer memory, the MLB approach can provide 
a complete and exact solution of the sign problem. We 
believe that similar ideas may also lead to the resolution 
of the sign problem in other fcrmion or real-time QMC 
schemes. 
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FIG. 2. Energy per electron, En/N, in units of Tiujo, for 
lo/a = 2 (squares) and lo/a = 8 (diamonds). Statistical errors 
are smaller than the symbol size. Open circles are taken from 
Ref. @ for lo/a = 2. 
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To conclude, we have proposed a multilevel block- 
ing approach to the fermion sign problem in finite- 
temperature QMC simulations. As presented, the 
method applies to the primitive path integral character- 
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